Ab initio complex band structure of conjugated polymers: 
Effects of hydrid DFT and GW schemes 
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The non-resonant tunneling regime for charge transfer across nanojunctions is critically dependent 
on the so-called (3 parameter, governing the exponential decay of the current as the length of 
the junction increases. For periodic materials, this parameter can be theoretically evaluated by 
computing the complex band structure (CBS) - or evanescent states - of the material forming the 
tunneling junction. In this work we present the calculation of the CBS for organic polymers using a 
variety of computational schemes, including standard local, semilocal, and hybrid-exchange density 
functionals, and many-body perturbation theory within the GW approximation. We compare the 
description of localization and /3 parameters among the adopted methods and with experimental 
data. We show that local and semilocal density functionals systematically underestimate the /3 
parameter, while hybrid-exchange schemes partially correct for this discrepancy, resulting in a much 
better agreement with GW calculations and experiments. Self-consistency effects and self-energy 
representation issues of the GW corrections are discussed together with the use of Wannier functions 
to interpolate the electronic band-structure. 



PACS numbers: 

I. INTRODUCTION 

The fields of molecular electronics and charge trans- 
port through nanojunctions have been deeply investi- 
gated in the past fifteen yearsP^ At the experimen- 
tal level many different techniques have been devel- 
oped, including those based on break junctions, nanos- 
tructured a nd scanning probe layouts, or self-assembled 
monolayers. 3 ' 4 Significant improvements in the accuracy 
with which these junctions are characterized have been 
achieved over the years, eg to address the I-V charac- 
teristics of single molecular junctions. Moreover, a large 
experimental literature exists^^ on non-resonant tunnel- 
ing experiments, where it is possible to determine the 
exponential decay (/?o) of the current I = Iq exp(—/3oL) 
as a function of the length L of the tunneling layer [such 
as eg a layer of organic self-assembled monolayer (SAM) 
connected to metallic electrodes] . Even though the /?o pa- 
rameter depends^J^ also on the detailed nature of the 
interface, it carries mostly information about the prop- 
erties of the tunneling layer itself, which makes /?o an 
interesting analysis and characterization tool. 

Experimentally, these measurements are performed us- 
ing different setups, ranging from metal- insulator- metal 
(MIM) junctions, as mentioned above, to the evaluation 
of kinetic constants of electro-transfer reactions (optically 



or electrochemically induced) in donor-bridge-acceptor 
molecular complexesP^H^ In terms of systems, measure- 
ments have been performed on a nu mber o f cases rang- 
ing from saturated olephins (alkanes P * 12 * 13 l to biological 
molecules (such as DNAG^HlIi) Theoretically, the elec- 
tronic mechanism underl ying the se experiments has been 
analyzed and understood! 1 * 5 * 8 * 1 -^ As stated in Ref. [5], the 
key parameter /? can be expressed (eg in MIM junctions) 
in terms of (i) the band gap E g and the (frontier) band 
widths (or hopping parameter) t of the insulating layer, 
and (ii) the alignment of the Fermi level in the metals 
with the energy gap of the insulator. Indeed, the effect 
of the electronic structure of the insulating layer can be 
singled out by evaluating^ the complex band structure 
(CBS), or evanescent states, in the limit of an infinitely 
long insulating region. The CBS approach is also particu- 
larly interesting for an ab initio evaluation of /?, where the 
calculations can be performed either using wavefunction- 
^ or Green's function-basecP approaches. A cartoon de- 
scribing the relation between the electronic structure of 
the MIM junction and the evaluation of the /3-decay fac- 
tor is given in Fig. [T] 

Nevertheless, since the f) parameter can be^ directly re- 
lated to the ratio between the energy gap and the band 
width of the insulator layer, the accuracy of standard 
electronic-structure simulations based on the Kohn-Sham 
(KS) framework of the density functional theory (DFT) 
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FIG. 1: The non-resonant tunneling experiment, (a) Scheme 
of the alignment of electronic levels in a metal-insulator-metal 
(MIM) junction. The complex and real band-structure (CBS 
and RBS) corresponding to the (extended) insulator system 
are reported in (b) and (c) respectively. The computed value 
to be compared with experiments is highlighted as /3(Ef), Ef 
being the Fermi energy of the MIM junction. 



{para phen ylene-vinylene) (PPV^PSland poly-(phenylene- 
imine}^2SI (PPI) , which are relevant from a technological 
point of view, where we compare also with recent ex- 
perimental dataP All chains are studied as isolated, see 
App. [Blfor full numerical details. 

II. METHOD: TRANSPORT 

To simulate the decay factors of non-resonant tunnel- 
ing experiments we adopt the CBS algorithm proposed 
by Tomfohr and Sankey (TS)P For a recent discussion 
of the connection of transport properties with the CBS 
theory see also the work by Prodan and CarP3 Within 
the TS approach, we need to evaluate the CBS in the 
limit of an infinitely thick insulating region (/? is in fact 
an asymptotic behavior). The outcome of this procedure 
is a set of (3(E) curves. The value /?o which has to be 
compared with the experiments is the smallest one (i.e. 
the most penetrating) aligned with the Fermi level of the 
junction: 



can be questioned. In fact, using eigenvalues computed 
from KS-DFT it is well known that the fundamental band 
gap is badly underestimated and (when using local and 
semilocal approximations) the delocalization of electronic 
states is typically overestimated. Moreover, the descrip- 
tion of this class of experiments in terms of single particle 
energies would require them to be interpreted as quasi- 
particle energies, in order to address the electronic dy- 
namics of the system. This is valid for advanc ed M BPT 
methods pSl such as Hedin's GW approximation}^^] and, 
at least in a perturbative sense, also for Hartree-Fock 
calculations. However, DFT Kohn-Sham states are ficti- 
tious orbitals with no direct physical interpretation, and 
their use in this context can only be justified by the as- 
sumption that the exchange-correlation kernel is an ap- 
proximation for the quasi-particle Hamiltonian. Model 
self-energies have also been recently^ used to correct 
the electronic structure of metal-molecule-metal junc- 
tions and found to be important for the evaluation of 
p. To date there has been no systematic investigation 
of the performance of different ah initio schemes in the 
calculation of the /3 decay factors. 

In this paper we study the effects of hy brid exchange- 
correlation functionals^ and the G W^ 1 * 22 ! approximation 
on the calculation of the /3 decay-factors, according to a 
scheme based on the the complex band structure (CBS) 
formalismP A detailed comparison with local and semilo- 
cal functionals is also provided. The theoretical back- 
ground of CBS, and GW and hybrid functionals are de- 
scribed in Sec. [IT|and |III| respectively. We discuss a gen- 
eral application of the Wannier functions interpolation 
to the case of GW electronic structure (Sec. IIIB). Our 
approach is applied to a number of polymers as reported 
in Fig. 2j We compute the CBS for poly-ethylene (PE) 
and poly-acetylene (PA) as references for saturated and 
conjugated chains respectively. We then consider poly- 



A) = P(Ef) 



(1) 



Since the electrodes are not considered in the calcula- 
tion, together with the proper metal-insulator interface, 
Ep is not known a priori and must be either estimated 
or calculated separately. This issue is discussed in de- 
tail in Ref. 8J. In the present work we will not com- 
pute the Fermi level alignment explicitly, and will rely 
on the estimation method proposed in the above refer- 
ence,^ which consists in evaluating (3(E) at the energy 



(a) (b) "f^K 




FIG. 2: Polymers studied: (a) PE (poly-ethylene), with 
each Carbon atom in the polymer chain being fully saturated 
by two Hydrogen atoms; (b) PA (poly-acetylene), where sp 2 
hybridization of Carbon atoms in the chain (each one also 
bonded to one Hydrogen atom) implies 7r-conjugation, i.e. 
alternation of single and double C-C bonds along the chain; 
(c) PPV [poly-(para phenylene-vinylene)], constituted by ben- 
zene rings connected through vinyl groups, also presenting 
conjugation; (d) PPI [poly-(phenylene-imine)], differing from 
PPV for the substitution of one Carbon atom in each vinyl 
group by a N atom, which still preserves the polymer conju- 
gation. 
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where df3/dE = (branch point). This method is orig- 
inally due to TersofP^ and based on the pinning of Ep 
by metal induced gap states (MIGS). We discuss this 
approximation in Sec. |V| where we compare with exper- 
imental results. Note that the value of j3 at the branch 
point is also connected (for non-metallic ID systems with 
a local potential) with the degree of lo calization of the 
density matrix, since /3 determine d 28 * 29 ! its spatial decay. 
Pictorially, a better description of 13(E) means then an 
improved description of the electronic localization. 

Scanning the energy spectrum, the CBS procedure 
searches for evanescent solutions to a given effective 
single-particle Hamiltonian. By definition, states with 
(real) Bloch symmetry k satisfy the relation 

f(R)Vk(r)=^ k (r + R) = AVk(r), (2) 

where R is any direct lattice vector, T(R) a translation 
operator, and A = e lk R . In the same way it is possible 
to define a complex Bloch symmetry k = k + i/3/2 setting 
A = e lK ' R = e lk R e~P' R / 2 , which is thus no longer a pure 
phase. The imaginary part of k implies a real space ex- 
ponential decay of the wavefunctions and it is customary 
to define (3(E) = 2 |Im[/«(i?) • e]| (e transport direction). 
The energy dependence of k comes from the fact that for 
a fixed energy E, the solutions are searched in terms of 
k, as it is usually done in scattering theory. 

By adopting a localized basis set {I&r)} (i and R 
orbital and lattice indexes, respectively), it is possible 
to define the Hamiltonian H and overlap S operators 
through their matrix elements ifjj(R) = {<f>io\H\(f)j R ) 
and Sij(H) = (<fiio\4>jti) , and the wavefunctions as \ipi) = 
EtR^ijO^) l&H.)- Setting Z = H — E S, the eigenvalue 
equation for H can be written as: 

JV 

Z(H m )C(R m ) = 0, (3) 

m=—N 

where matrix multiplication is implicit, and N is the 
number defining the last non-zero matrix Z(Rjy). Here 
we are assuming a real space decay of the Hamiltonian 
and overlap matrices, which is typically physical even if 
the range can be strongly dependent on the scheme used 
to define the Hamiltonian. We will comment later on this 
point when discussing the use of GW or HF methods. It 
is then possible to derive^ the following system of 2N 
matrix equations: 

N-l 

- Z(R m )C(R m ) = Z(R N )\C(R N - 1 ), (4) 

m=-N 

C(H m+1 ) = AC(R m ), m = —N, N — 2. (5) 

where Eq. ^ is a direct consequence of Eq. Q. Such 
matrix equations become an eigenvalue problem for A as- 
suming we can invert the matrix Z(Rn). As pointed out 
in Ref . |8 , this is very often not the case as the matrix is 
singular, but the singluarities can be avoided. The reader 



is referred to the original work for the details. The full al- 
gorithm pro posed in Ref. [8] has been implemented in the 
WanT cod e 1 30 * 31 -! and used for the present work. A sim- 
ple tight-binding analytical model (a generalized version 
of the one presented in Ref. [5]) is discussed in App. [A] 
This model will be used in Sec. |IV| to fit and interpret 
the real and complex band-structures of the polymers we 
have investigated here. 



III. METHOD: ELECTRONIC STRUCTURE 

According to the above discussion, in order to simulate 
the decay coefficient of a MIM junction, we need to com- 
pute the electronic structure of the insulating layer (con- 
sidered as infinitely extended). The underlying reason 
for this simplification is that we interpret the computed 
single-particle energies of the system as quasi-particle 
(QP) energies, which in turn determine the dynamics 
of singly charged excitations. In general the transport 
problem for interacting systems is more complicated than 
that and requires a more sophisticated treatment 
For DFT, the understanding of how the exact KS-DFT 
Hamiltonian performs to compute transport has been re- 
cently the subject of several invest igat ions BUSH Apart 
from the properties of the exact functional, currently 
available DFT approximations like LDA or GGA have 
been demonstrated to systematically overestim ate the 
conductance, especially for off-resonant junctions! 32 * 45 * 46 ^ 
Pragmatically, this suggests that corrections beyond lo- 
cal and semilocal KS-DFT approaches are needed. 

Using QP-corrected electronic structure to compute 
charge transport (eg by means of the Landauer for- 
mula 4 -^) through interacting systems seems to be a rea- 
sonable approximation when finite-lifetime effects are 
weakP2E3 Indeed, a number of works com puting Q P 
energies by mean s of t he GW approximatiorP^E3EU or 
model self-energies^E2HHl have been reported in the lit- 
erature. On the other hand, hybrid-exchange functional 
methods like B3LYP^ or PBEG^S are widely used and 
lead to band gaps and band widths which are usually 
closeJ^to the experimental values than simple semi- local 
KS approaches, for both molecules and solids. Recent 
workgEiHSZ] have further investigated the accuracy of hy- 
brid exchange functionals, also in comparison with GW 
calculations. In this work we compare GW and hybrid 
functionals for the calculation of electronic structure and 
transport properties of selected organic polymers. In the 
following we summarize the GW approximation and un- 
derline some formal similarities with hybrid functionals. 



A. The GW approximation and hybrid DFT 

A many-body theoretical formulation of the electronic 
structure problem can be obtained by using the Green's 
function formalism. The one-particle excitation energies 
of an interacting system are the poles of its interacting 
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Green's function G(£))22IS2I w hi c h can be written as: 

G(E) = [EI - h Q - £(£)]- 1 , (6) 

where ho is an effective single particle Hamiltonian and 
is the non-local, non-hcrmitcan, frequency depen- 
dent self-energy operator. In general, X is not known a 
priori and must be approximated. In this work th e self- 
energy is computed within the GW approximation's^ 

Z GW (r 1 ,r 2 ,E) = i J ^e-^'Gin^E-u 1 ) x 

W^n.ra.w'), (7) 

where W(lu) is the screened Coulomb interaction eval- 
uated in the random phase approximation (RPA). For 
more details see e.g. Refs. |63|64j . In the simplest im- 
plementation of the GW approximation, the self-energy 
is computed non self-consistcntly, i.e. by evaluating G 
and W according to the eigenvalues and eigenvectors of a 
reference non-interacting Hamiltonian (typically the KS 
Hamilotonian at the LDA or GGA level). Such a proce- 
dure is known as Go Wo and gives reasonable results f or 
the quasiparticle energies in a number of cases ESESESI i n 
the present paper we exploit the Go Wo approximation, 
and evaluate the frequency integrals in Eq. Q by using 
a plasmon pole model according to Godby and NeedsP^l 
In this work, the main quantities we are interested in 
are the QP energies. If we neglect finite lifetime effects 
and take (or symmetrize) the self-energy to be hermitean, 
QP energies are given by first order perturbation theory 
as: 



KS 



+ (i> m \Z(e% p )-v xc \^ m ), (8) 



As customary,^ in order to solve for e^ p , the self-energy 
in the above equation is expanded to first order as a func- 
tion of E. 

In order to get mo re ph ysical insight we also refer to 
the (static) COHSEXPS] approximation of GW, where 
the self-energy is written as £ = £ COH + £ SEX : 

£ SEX (r 1; r 2 ) = -Ttri.rajW^ra.O), (9) 
E COH (n,r 2 ) = ^d(r 1 ,r 2 )W p (v ll r 2 ,0). (10) 

Here 7(1*1, r 2 ) is the one particle density matrix, and 
W p = W — v is the dynamical contribution to W (be- 
ing v the bare Coulomb interaction). The COHSEX self- 
energy is thus the sum of a statically screened exchange 
term and a local potential. This partition is particu- 
larly useful for discussing the connection between hybrid 
exchange functionals and GW. In the former case, the 
potential can be written asP^l 



,hyb 



a v 



NL 



(11) 



where is the non-local exchange potential, while 
and are local potentials. It is then straightforward to 
interpret a as an inverse effective screening to stress the 



formal analogy of Eqs. ([9-10 1 and (11). Similar consid- 
erations are of course valid for more complex forms of 
hybrid functio nals, like range separated or local formu- 
lations .HmUMl This formal analogy is well-known in the 
literature,^ and it has also been further investigated re- 
cently.^ Besides their accuracy for thermochemistry, this 
analysis highlights that also the electronic structure com- 
puted by non-local hybrids can benefit from the inclu- 
sion of some screened exchange term. Indeed, improved 
description of the electronic struct ure f or finite and ex- 
tended systems are typically foundp^^ even though the 
accuracy may vary significantly depending on the system. 



B. Interpolation of GW electronic structure using 
Wannier functions 



Dealing with periodic systems, interpolation over the 
first Brillouin zone (1BZ) is a long standing issue. Cal- 
culations are typically performed by discretizing k-points 
in the 1BZ, and some post-processing schemes [such as 
eg those to compute density-of-states (DOS), Fermi sur- 
face, band structure, or phonons] might need a better 
discretization of 1BZ than some of the previous steps (of- 
ten aimed at computing total energy, forces, and charge 
density). This is particularly critical when the compu- 
tational requirements of the adopted methods limit the 
k-point discretization, as is the case for GW calculations. 
Schemes able to refine or interpolate^ over the 1BZ are 
particularly useful for this pur pose. One of these is the 
Wannier interpolation? ^ 2 "^ where the localization of 
the Wannier function (WF) basis together with the fi- 
nite range in real space of the Hamiltonian are used to 
perform a Fourier interpolation of the eigenvalues, and 
eventually eigenvectors. The use of this scheme to inter- 
polate GW results has been also reported elsewhere by 
Hamann and Vanderbilt.^ 

The procedure can be applied not only to the Hamil- 
tonian, but in principle to any operator A(r, r') with the 
translational symmetry of the Hamiltonian. First, we 
define the projector P over a subspace of interest: 



P = J2 IV-nO^kl 



(12) 



nk 



in terms of the eigenvectors of H . When the subspace P 
is complete, we can represent A as: 

A = X^|^ m kM mn (k)(V>„k| (13) 
k mn 

A mn (k) = (^ mk |^|^„k) (14) 

Here, A is diagonal with respect to the k-index because 
it commutes with the translation operators of the direct 
lattice (as assumed). In practice, limiting the number 
of eigenstates of H included in P is equivalent to con- 
sidering the projection of A on the P subspace, namely 
A p = PAP instead of A. At this point we can use 
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the definition of maximally localized Wannier functions 
(MLWF's), according to Ref. [73] : 



(15) 



to obtain an expression for the matrix elements of A on 
the Wannier basis, AfAH) = (wio\A\wjR): 



(16) 



Note that when the original Marzari-Vanderbilt proce- 
dure^ is applied without any disentanglement the U k 
matrices are a unitary mapping of N Bloch states (usu- 
ally occupied, but not necessarily) into N WF's. Instead, 
when the disentanglement is performed, the resulting WF 
set does not span the whole P subspace (t/ k are then rect- 
angular matrices). This means that in the general case 
the final representation of A is actually not projected on 
P but on the smaller subspace spanned by the WFs. 

Assuming that A p is decaying fast enough in real space 
to have ||A P (R)|| ~ for |R| > |R | (where R is within 
the finite set compatible with the initial k-point grid) , we 
can perform the following Fourier interpolation to obtain 
the matrix elements of A for any k' point: 



A^(k') = 



|R|<|Ro| 

E 

R 



ik'R 



Ag(R) 



(17) 



When interpolating GW results, we want to represent 
the operator £(P) = Y, GW {E) — w xc which is in general 
non-local, non-hcrmitean and frequency dependent. For 
the sake of the Wannier interpolation, we are mainly in- 
terested to check that the intrinsic non-locality of PEP 
is compatible with the selected k-point grid (or, in other 
terms, that the GW calculation is converged with respect 
to the number of k-points used) . The localization of the 
GW self-energy is further discussed in Sec. |V A| espe- 
cially in connection with the usual approximation that 
neglects off-diagonal E m „(kP) matrix elements. 



C. Numerical approach 

In this work, DFT and hybrid-DFT calculations have 
been performed using the CRYSTAL09 packageP The 
code implements all-electron electronic structure meth- 
ods within periodic boundary conditions and adopts an 
atomic basis set expanded in Gaussian functions (further 
details in App. |B| , Once the Hamiltonian matrix ele- 
ments are obtainedp' the real and complex band struc- 
tures are interpolated (as di scussed in the previous Sec- 
tions) using the WanT3M3 package. 

GW results have been obtained using the pla ne- wave 
and pseudopotentials implementation of SaXP^I which is 
interfaced to Quantum-ESPRESSC 2 ^ (QE) for DFT 
calculations. In this case, once the Kohn-Sham electronic 



TABLE I: Lattice parameter c [A] for PA, PE, PPV, and 
PPI, computed using different XC schemes as implemented 
in CRYSTAL09. In the case of PA, we also report the bond 
length alternation [A] (BLA) of single and double C-C bonds. 



Scheme 


PA 


PA-BLA 


PE 


PPV 


PPI 


LDA 


2.463 


1.369/1.416 


2.537 


6.644 


12.869 


PW 


2.482 


1.376/1.428 


2.570 


6.712 


13.001 


BLYP 


2.493 


1.380/1.435 


2.590 


6.747 


13.079 


PBE 


2.484 


1.378/1.429 


2.572 


6.719 


13.014 


B3PW 


2.471 


1.362/1.432 


2.560 


6.692 


12.961 


B3LYP 


2.476 


1.375/1.435 


2.570 


6.706 


12.998 


PBEO 


2.468 


1.359/1.432 


2.553 


6.680 


12.938 


HF 


2.465 


1.332/1.457 


2.556 


6.689 


12.950 



structure is evaluated, we first compute MLWF' j 73 * 74 ! us- 
ing WanT and then apply the CBS technique. In or- 
der to assess any systematic error in comparing GW and 
hybrid-DFT results (which have been obtained using dif- 
ferent basis sets such as plane waves and local orbitals), 
we have also performed hybrid-DFT and HF calculations 
using QE and SaX. In this case WFs are computed on 
top of the already corrected electronic structure. Results 



are shown for the case of polyacetylene (see Tab. Ill and 
Fig.[8^b) in particular). The excellent agreement between 
the two sets of data suggests that the pseudopotential ap- 
proximation, the basis set, and the numerical thresholds 
are sufficiently well converged to have negligible influence 
on the results presented. Full computational details and 
parameters are reported in App 15] 



IV. RESULTS 
A. Structural properties 

Before focussing on the electronic and transport prop- 
erties of the isolated polymer chains from Fig. [2j we 
investigate their structure by fully relaxing both the 
atomic positions and the cell parameters using different 
exchange-correlation schemes. All systems are treated 
with one-dimensional periodicity, the details of the calcu- 
lations (performed using CRYSTAL09) have been given 
in App. [B] The results for the lattice parameters are re- 
ported in Tab. [T] 

In the case of PA, electronic properties such as the 
band gap (as well as the evanescent states) are strongly 
dependent on the dimerization of the C-C bond lengths 
(Peierls distortion). Such bond alternation is not eas- 
ily captured by local and semilocal DFT schemes leading 
to band gaps that are far too small, i.e. to an overes- 
timation of the metallicity. Since these parameters are 
critical to our purpose, we report also the bond length 
alternation (BLA) for PA in Tab. [!] Our HF results 
are in excellent agreement with previously published re- 
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suits PflSU l n the following we will label as PAhf the cal- 
culations performed using the geometry from Ref. |80j . 
where c = 2.469A and BLA is 1.339/1.451 A. We also 
decided to consider two frozen geometries of PA (namely 
PAi and PA2) according to other data in the litera- 
turePsHln the case of PAi^we set c = 2.451 A and BLA 
to 1.370/1.460 A, while for PA^c = 2.496 A and BLA to 
1.340/1.540 A. In passing we note that the PAi geometry 
is also very similar to the one adopted in Refs. [8"3"H8"S] . 
where c = 2.451 A and BLA is set to 1.360/1.440 A, ac- 
cording to experimental data!££123 Since these theoretical 
studies report GW results, in Sec. |IVB] we will compare 
with PAi. 

Data from multiple geometries of PA are useful to de- 
couple the electronic and structural effects of the adopted 
XC schemes on the complex band structure. Since the 
extent of this goes beyond PA, in the following we de- 
cided to look at the effect of XC treatment first using a 
frozen geometry, independently of the adopted method, 
and then to compare also with the same quantities ob- 
tained using fully relaxed geometries. Moreover, since 
GW corrections are usually computed without perform- 
ing a further structural relaxation, working at fixed ge- 
ometry allows us to compare GW corrections and results 
from hybrid functionals for identical geometries. For each 
polymer except PA we adopted the geometry obtained 
after full relaxation at the PBE level using Quantum- 
ESPRESSO. Lattice parameters (2.564 A for PE, 6.702 
A for PPV and 13.004 A for PPI) are in very good agree- 
ment with those in Tab.|T]obtained using PBE in CRYS- 
TAL09. 



B. Electronic structure: poly-ethylene and 
poly-acetylene 

In this section we investigate the effect of different XC 
schemes (local and semilocal DFT, hybrid functionals, 
HF, and GW) on the electronic properties of two proto- 
type polymers (PE and PA), while results for PPV and 
PPI will be reported in Sec. |IV C| We compute both the 
real and the complex band structures. Figures[3]|4]refer to 
the case of frozen geometries (i.e. geometry is not chang- 
ing according to the adopted scheme, see also Sec. IV A I. 
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Details about the electronic structure of the three PA 
geometries studied are reported in Tab. |III| In the case 
of PA (as well as PPV), we can also comp are with pre- 
viously published GW calculationsPHS51S2l Note that the 
GW results are interpolated using MLWF's, which re- 
sults in filtering out some of the states above the vacuum 
level. 

In Fig. [3], we reportP^the real and complex band struc- 
ture for PE. Our results are in reasonably g ood agreem ent 
with previously published theoretical dataP^* 19 * 89 * 90 ! As 
can be seen from panels (a,b), the maximum value of 
from the arc across the fundamental gap is not chang- 
ing much when passing from PBE to PBE0 (/3 max — 0.8 
A -1 ). In agreement with Ref. [TU], we believe this to be 




FIG. 3: Poly-ethylene (PE). Real and complex band struc- 
ture using the schemes: (a) PBE, (b) PBE0, (c) HF and (d) 
GW. Solid (black) lines in (a,b,c) refer to CRYSTAL09 cal- 
culations. GW results are obtained using SaX (circles), and 
interpolated with WFs (solid black lines). 



one of the reasons why the /? computed at the LDA/GGA 
levels have been found in agreement with the experiment. 

The GoWo-corrected band structure is reported in 
panel (d). It compares favorably with existing literature 
dataP^The main qualitative difference with Fig. 2(a) of 
Ref. [19j comes from the filtering of some vacuum-like 
states, due to the use of the localized basis (or WF inter- 
polation) in our approach. Our results are more similar 
to those presented in Ref. [8], obtained using a localized 
basis set. This has little influence on the part of the CBS 
spectrum that is physically relevant for the non-resonant 
tunneling experiment. The GqWo CBS is not computed 
here because the missing self-consistency is critical for 



CBS, as it is discussed in Sec. V A Instead, we have com- 
puted CBS from the self-consistenPUCOHSEX electronic 
structure. All results for /3 max are collected in Tab. [TTJ 
where we compare data obtained by using CRYSTAL 



TABLE II: PE: maximum of /3(E) [A -1 ] inside the gap 
computed by means of different theoretical schemes. Re- 
sults are obtained by CRYSTAL09 (CRY) and by Quantum- 
ESPRESSO (QE)-SaX. 



/^max 






poly ethylene 


(PE) 






LDA 


PBE 


PBE0 


HF 


COHSEX 


CRY 


0.77 


0.77 


0.81 


0.92 




QE-SaX 


0.80 


0.81 
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0.94 
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FIG. 4: (color online). Poly-acetylene (geometry PAi). Real 
and complex band structure using the schemes: (a) PBE, (b) 
PBE0, (c) HF and (d) GW. Solid (black) lines in (a,b,c) refer 
to CRYSTAL09 calculations. GW results are obtained using 
SaX (circles), and interpolated with WFs (solid black lines). 
Dashed (red) lines refer to data interpolated according to the 
tight-binding model in App. [A"] 



and QUANTUM-ESPRESSO or SaX after Wannieriza- 
tion. Results and trends compare reasonably well. 

In Fig. |4j we show the results for real and complex 
band structure in the case of PAj . The dependence of the 
ft parameter on the exchnage-correlation scheme is dis- 
played, exchange and correlation for the j3 parameter is 
displayed. As expected, when increasing the percentage 
of non-local exchange from PBE (0%), to PBE0 (25%) 
up to HF (100%), the band gap opens, and the "degree 
of localization" increases as indicated by the increasing 
values of /3 max . This trend is general and also found for 
the other conjugated polymers studied in this work. A 
detailed description of the computed band structure for 
PA in the PAhf 7 PAi, and PA2 geometries is reported in 
Tab.|nT]for all the methods. 



In order to gain more physical insight from our calcu- 
lations, we have also fitted the data by using the gener- 
alized nearest-neighbors (NN) model presented in Sec. [TT] 
and App. [Xj This model has three parameters, which 
approximately correspond to the band gap E g , and the 
band widths of the HOMO and LUMO bands (related to 
the parameters t\ and £2)- The proper relation between 
the band widths and t\ 2 is given in Eqs. (All -A12|. 
The main difference between this model and the one in- 
troduced in Ref. [8] is that the one used here allows for 
different band widths for the HOMO and LUMO. While 
this difference is found to be small (but generally not neg- 
ligible) in the cases studied, the generalized model allows 



TABLE III: PAhf, PAi, PA 2 : Band gap (E g , [eV]), hop- 
ping paramete rs (t i, t2, [eV]) according to tight-binding (TB) 
model in Eq. |A5|, maximum of /3(E) [A -1 ] inside the gap. 



Calculations are performed using DFT, hybrid-DFT, HF, and 
diagonal-GW schemes. /3 max through the CBS interpolation 
by using the TB model are also shown. All the data have 
been computed using CRYSTAL09 except the GW results 
and the X-QE lines (X=LDA,PBE,PBE0). 



Scheme 


Eg 


ti 


t 2 


/3 m ax 


/omodcl 
Pmax 


poly-acetylene (PAhf) 


LDA 


0.99 


2.97 


0.199 


0.13 


0.13 


PW 


1.01 


2.97 


0.203 


0.14 


0.14 


BLYP 


1.02 


2.96 


0.199 


0.14 


0.14 


PBE 


1.01 


2.98 


0.204 


0.14 


0.14 


B3PW 


1.94 


3.78 


0.223 


0.19 


0.21 


B3LYP 


1.95 


3.77 


0.218 


0.20 


0.21 


PBE0 


2.20 


3.98 


0.228 


0.21 


0.22 


HF 


6.97 


6.45 


0.277 


0.38 


0.43 


d-GoWo 


2.60 


3.85 


0.154 


(0.17) 


0.27 


LDA-QE 


0.99 


3.02 


0.214 


0.13 


0.13 


PBE-QE 


1.02 


3.00 


0.212 


0.14 


0.14 


PBE0-QE 


2.22 


3.99 


0.238 


0.20 


0.22 






poly-acetylene (PA 


1) 




LDA 


0.78 


2.90 


0.144 


0.11 


0.11 


PW 


0.80 


2.90 


0.148 


0.11 


0.11 


BLYP 


0.80 


2.89 


0.144 


0.11 


0.11 


PBE 


0.80 


2.91 


0.149 


0.11 


0.11 


B3PW 


1.64 


3.73 


0.166 


0.17 


0.18 


B3LYP 


1.64 


3.72 


0.161 


0.17 


0.18 


PBE0 


1.88 


3.94 


0.171 


0.18 


0.19 


HF 


6.46 


6.46 


0.212 


0.36 


0.40 


d-GoWo 


2.05 


3.72 


0.090 


(0.14) 


0.22 






poly-acetylene (PA 


2) 




LDA 


1.68 


2.76 


0.134 


0.24 


0.24 


PW 


1.72 


2.76 


0.139 


0.25 


0.25 


BLYP 


1.72 


2.75 


0.134 


0.25 


0.25 


PBE 


1.72 


2.77 


0.139 


0.25 


0.25 


B3PW 


2.89 


3.49 


0.153 


0.31 


0.33 


B3LYP 


2.90 


3.47 


0.149 


0.31 


0.33 


PBE0 


3.20 


3.67 


0.158 


0.33 


0.35 


HF 


8.37 


5.88 


0.196 


0.50 


0.56 


d-GoWo 


4.17 


3.87 


0.092 


(0.27) 


0.43 



for a more accurate fitting of the electronic structure of 
polymers. As can be seen in Fig. Qa-c), the model fit- 
ting (dashed, red lines) is very accurate for all the local, 
semilocal and hybrid functionals. In general the HF data 
shows the largest deviation from the model. We believe 
the non-local nature of the exchange potential to be the 
origin of this behavior. 

Figure |4|d) reports the GW results for PAi. The 
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FIG. 5: (color online). PPV: Real and complex band struc- 
ture, as in Fig. [4] 



fundamental gap is 0.8 eV at the PBE KS-DFT level, 
while it increases to 2.05 eV when GqWq is applied. 
This is in very good agreement with previous GW data 
for pApHUMl As already mentioned, GW calculations 
are performed by using plane-waves and pseudopoten- 
tials, while hybrid-DFT is evaluated on a localized basis 
set. In order to assess a possible systematic error due 
to this procedure, we have also computed the electronic 
structure of PA H f at the LDA, PBE and PBE0 levels, 
by using the plane-wave implementation of Quantum- 
ESPRESSO(QE). Results are reported in Tab. [TTT] at 
the lines X-QE (X=LDA,PBE,PBEO). 



The two sets of 
data are found to be in excellent agreement, allowing for 
a direct comparison of GW and hybrid-DFT data. See 
also Sec. VA and Fig. [8lb) for further details. 



As a last remark, acco rding to the standard approach 
for GW calculations j 63 * 64 ! only the Kohn-Sham eigenval- 
ues are corrected, without modifying the wave functions. 
In other words, only the diagonal matrix-elements of the 
self-energy (on the original DFT Bloch eigenvectors) are 
considered. The effects of this is analyzed in detail in 



Sec. V A when it will be demonstrated that this proce- 
dure has a sizable effect on the calculation of the complex 
band structure. For this reason, the GW-CBS directly 
computed is reported in a lighter color in the right panel 
of Fig. Qd) , while /3 max is put in parentheses in Tab. Ill 



C. Electronic structure: PPV and PPI 

In this Section we consider two further polymers, 
namely poly-para phenylene-vinylene (PPV) and poly 
phenylene-imine (PPI). PPV has been largely investi- 



TABLE IV: PPV: Band gap (E g , [eV]), hopping parameters 
(ii, t2, [eV]) maximum of /3(E) [A -1 ] inside the gap (com- 
puted and fitted from the tight-binding model), as in Tab. Ill 
Last column reports the value /3 max for the relaxed geometries 



Scheme 


Eg 


ti 


t 2 


/3max 


/□model 
Pmax 


/□relax 
rmax 




poly para-phenylene 


-vinylene (PPV) 




LDA 


1.28 


1.07 


0.016 


0.19 


0.18 


0.18 


PW 


1.31 


1.07 


0.017 


0.19 


0.18 


0.19 


BLYP 


1.31 


1.07 


0.016 


0.19 


0.18 


0.19 


PBE 


1.31 


1.07 


0.017 


0.19 


0.18 


0.19 


B3PW 


2.22 


1.39 


0.021 


0.25 


0.23 


0.26 


B3LYP 


2.22 


1.38 


0.020 


0.25 


0.23 


0.26 


PBE0 


2.46 


1.46 


0.022 


0.26 


0.24 


0.28 


HF 


6.74 


2.46 


0.033 


0.41 


0.38 


0.46 


G Wo 


3.09 


1.63 


-0.004 


0.20(d) 


0.28 
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FIG. 6: (color online). PPI: Real and complex band struc- 
ture, as in Fig. [4] The fitting procedure has been applied on a 
cell with half length and then folded, resulting in four bands 
instead of two. 



gated for its role in organic (opto-)electronicspSI while 
oligo phenylene-imine molecules attached to gold leads 
have been recently considered and the j3 decay coefh- 
cents measured This makes these two polymers par- 
ticularly appealing for our analysis. 

Results for PPV are reported in Fig. [5] and Tab. [TV] 
The behavior of the real and complex band structures are 
qualitatively in agreement with what we have found for 
PA. PBE0 results (/3 max = 0.26 A" 1 ) are those among 
the hybrid functionals that best compare with the in- 
terpolated GW data (/3 max = 0.28 A^ 1 ), though slightly 
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TABLE V: PPI: Band gap (E g , [eV]), hopping parameters 
(ti, ti, [eV]) maximum of P{E) [A -1 ] inside the gap (com- 
puted and fitted from the tight-binding model), as in Tab. Ill 
Last column reports the value /3 max for the relaxed geometries 



Scheme 


Eg 


ti 


ta 


/3 m ax 


/□model 
Pmax 


/□relax 
Pmax 


poly-phenylene-imine (PPI) 


LDA 


1.40 


1.05 


0.014 


0.21 


0.20 


0.20 


PW 


1.42 


1.05 


0.014 


0.21 


0.20 


0.21 


BLYP 


1.42 


1.04 


0.014 


0.21 


0.21 


0.21 


PBE 


1.42 


1.05 


0.014 


0.21 


0.20 


0.21 


B3PW 


2.38 


1.36 


0.020 


0.27 


0.26 


0.29 


B3LYP 


2.38 


1.35 


0.019 


0.27 


0.26 


0.29 


PBEO 


2.64 


1.43 


0.021 


0.29 


0.28 


0.30 


HF 


7.12 


2.38 


0.036 


0.45 


0.43 


0.50 



underestimating /3 max and to a larger extent the band 
gap. In the case of PPI, results are reported in Fig. [6] 
and Tab. |Vj For this polymer, the model fitting has been 
applied on a cell with half length and then bands have 
been folded leading to four interpolated bands instead 
of two. Despite the reduction of translational symme- 
try, this procedure leads to a better fit because it de- 
scribes a larger part of the frontier electronic structure. 
In agreement with previous cases, while we do not have 
GW results for PPI, we consider PBEO data (/3™^ x =0.29 
A -1 ) as our best estimate. In the next Section we will 
also discuss the comparison of th ese computed data with 
recent experimental results.^H For the PPV and PPI 
cases, we have also studied separately the effect of the 
geometrical relaxation induced by the different function- 
als on the CBS. Such effect is consistent with the trends 
already observed at fixed geometry. In the last column of 
Tabs. |IV[[V| we report the /3 max value for the relaxed poly- 
mer geometries (labelled as /^ax*)- The /3 max parameters 
increase further with increasing fraction of non-local ex- 
change, when the geometries are relaxed according to the 
adopted functional. While the coupling of the electronic 
structure with the structural properties is particularly 
evident and critical in the case of PA (since the opening 
of the gap is due to Peierls distortion of the C-C bonds) , 
it is much less pronounced for PPV and PPI where it 
accounts for a correction term only, the leading contri- 
bution being the description of the electronic levels. 



V. DISCUSSION 

A. Analysis of the GW data 

In this Section we discuss the effect of some of the ap- 
proximations involved in the evaluation of the GW self- 
energy. In particular we address issues related to the rep- 
resentation of £ when computing the CBS, as well as the 



effect of the self-consistency. The GW self-consistency 
is investigated within the static COHSEX approxima- 
tion. In doing so, we discuss the localization properties 
of the resulting Hamiltonians together with the quality 
of the NN model fitting of the CBS. We focus on the 
case of PA, which is a good prototype for this class of 
one-dimensional systems. 

Let us begin with the representation problem. As al- 
ready recalled in Sec. 



IIIB 



assuming a large enough sub- 
set of Bloch vectors (in principles, all of them) , the self- 
energy operator can be represented as: 



(18) 



[see also Eqs. ( 13|14 |], In the usual GW practice, be- 
sides evaluating the self-energy by using the underlying 
KS-DFT electronic structure for G and W (GqWo ap- 
proximation, i.e. no self-consistency), it is also custom- 
ary to neglect the off-diagonal band indexes m / n in 
Eq. (181 when computing QP energies. This approxima- 



tion forces the self-energy to be diagonal on the KS-DFT 
Bloch states, and thus allows us to modify the quasi- 
particle energies without changing the DFT wavefunc- 
tions. If we assume that the correctly represented self- 
energy [Eq. ( 18 )] is physically short-ranged in real space 
(further comments follow), as it happens eg for HF and 
COHSEX, when the representation is taken to be diag- 
onal on the Bloch basis spurious long range components 
of the self-energy may (and typically do) arise, as evi- 
dent from simple Fourier transform arg uments. The di- 
agonal approximation has been found^ to have little ef- 
fect on the quasi-particle energy spectrum (at least when 
LDA wavefunctions are a reasonable starting point). Our 
investigations confirm this picture at both the HF and 
COHSEX level. As discussed in the following, the case 
of the complex band structure is more critical. 

First we focus on a tight-binding model. In Fig. [7] we 
report the real and complex band structures for such a 
model, according to App. [XJ In order to simulate the ef- 
fect of a diagonal self-energy, we refer to a picture where 
the £ correction can be modelled as a stretching of the 
bands (which may be different for valence and conduction 
states) plus a scissor operator applied to the HOMO- 
LUMO gap. A simple scissor and a scissor+stretching 
corrections are applied to the model in Fig. 7]^a,b) (solid 
black lines for the original model, thick light gray lines 
for the corrected ones). The new Hamiltonian includ- 
ing the corrections is now longer ranged than the orig- 
inal nearest neighbors Hamiltonian h°, because of the 
non-local projectors used to express the scissor and scis- 
sor+stretching corrections. The different spatial decay of 
the pristine and corrected Hamiltonians is shown in the 
lower panels of Fig. [7j We then fit the corrected Hamil- 
tonian by using again the NN model (dashed, red lines). 
This fitting emulates the real band structure, but using 
a short ranged (nearest-neighbors) Hamiltonian. The ef- 
fect on the CBS is evident and sizeable. The simply 
shifted and stretched electronic structures leads to little 



10 




(a) HF Q (LDA) 

0.2 0.4 



FIG. 7: (color online). Real and complex band structure for 
a nearest-neighbors (NN) tight binding model Hamiltonian 
h° (thin solid black line). A scissor and scissor+stretching 
corrections to h° are applied and shown in panels (a) and (b), 
respectively. The thick solid gray lines represent the elec- 
tronic structures obtained while including such corrections. 
The dashed red lines are obtained by a NN tight binding 
fitting of the real electronic structure after the corrections. 
Lower panels report the spatial decay of the original and cor- 
rected Hamiltonians. 



corrections to the CBS, while much bigger corrections 
are obtained when considering the fitted short-ranged 
Hamiltonians. Because CBS measures the decay of the 
evanescent states in real space, it is not surprising that a 
method which does not update the wavefunctions (as the 
diagonal self-energy corrections) is not able to capture 
all of the physics involved in the change of the electronic 
structure. 

In order to further numerically support this interpre- 
tation and to investigate the effect of the self-consistency, 
we have evaluated the HF and COHSEX self-energies 
for PAjjf- At first we have done so non-self-consistently 
(self-energies evaluated on the LDA wavefunctions) , with 
and without the diagonal approximation. Then self- 
consistentPH COHSEX results are provided. HF results 
are plotted in Fig. [8]Ja,b), COHSEX data in panel (c). 
LDA /3 max is shown with a dashed line in the CBS pan- 
els as a reference. Regarding the real band structure, we 
found that the inclusion of off-diagonal matrix elements 
is not very relevant for PA, the bands being almost over- 
lapping for both HF and COHSEX [diagonal data shown 
in panels (a,c) by dashed gray lines]. The situation is 
different for the CBS, as it is highlighted in Fig. ^a,). 
The HF correction of /3 max from the full off-diagonal 
representation is almost twice as big as the diagonal cor- 
rection. This confirms the behavior observed with the 
models in Fig. [7] 

This observation also correlates with the decay of the 
HFo(LDA) Hamiltonian reported in the lower part of 
panel (a). As for the models [Fig. [7], the diagonal rep- 
resentation induces a much longer (and unphysical) de- 
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FIG. 8: (color online). Real and complex band structure for 
PAhf at the HF and COHSEX levels. Panel (a) shows non- 
self-consistent HF results with a diagonal (gray lines) and 
fully off-diagonal (black lines) representation of the correc- 
tion. Panel (b) reports the same data for the self-consistent 
HF solution, as computed from SaX (solid black lines) and 
CRYSTAL (green triangles). Panels (c) shows COHSEX 
data: Diagonal non-self-consistent (gray lines) and fully self- 
consistent (black lines) data are reported. In both cases, a 
NN tight-binding fit of the real and complex band structures 
is performed (dashed red lines). Lower panels show a mea- 
sure of the spatial decay of the COHSEX and HF Hamiltonian 
matrices on the WF basis. 



cay. The same situation is found for COHSEX (the off- 
diagonal results at the first SCF iteration are not shown) . 
The proper Hamiltonian decay (black line, circles) is 
clearly longer ranged than the LDA results, because of 
the non-local contribution of the exchange operator. The 
decay of the exchange potential is driven by that of the 
density matrix, which in turn is related^ to (3 at the 
branch point. Being /3 max typically underestimated at 
the LDA level, so is the decay of the HFo Hamiltonian. 
The effect of self-consistency of HF (as well as COHSEX) 
is then to reduce such over-delocalization and to produce 
shorter ranged self-energies. This is shown in the decay 
plot of panels (b,c). 

This behaviour has strong consequencies regarding the 
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quality of the NN model fit. In the case of local and self- 
consistent hybrid functional calculations from CRYS- 
TAL, the model fit (dashed red lines in Figs.[4j[5| works 
very well when compared with the full calculations, de- 
spite its semplicity. This is not the case when comparing 
with the diagonal corrections of Fig. |8^a,c). The off- 
diagonal representation improves the situation but does 
not solve the problem. The failure of the model fit in 
Fig. [8](a) is in fact mostly related to the decay of the ex- 
change operator. A much better fit is then found when 
full self-consistency is included, as in Fig.[8pD,c). Finally, 
we have also compared the HF results from SaX (plane- 
waves and pseudopotentials) and CRYSTAL (localized 
basis, full electron). Results are reported in Fig.|8jb) by 
the black solid lines and the green triangles, respectively. 
We find excellent agreement for both the real and com- 
plex band structures, confirming that the results from the 
two codes are well comparable and free of any systematic 
error. 

In the light of the discussion above, when pre- 
senting the diagonal Go Wo data for PA and PPV 
[Figs.gd) jljd)], the CBS is better described by the inter- 
polated data (dashed lines) instead of that directly calcu- 
lated from the diagonal GW corrections (solid thin lines) . 
Similar conclusions about the importance of describing 
changes to wavefunctions when applying GW to trans- 
port calculations are reported also in Refs. |49I50I93| . 



B. Electronic structure 

Now we can turn to discussing the accuracy of the 
electronic structure calculations for conjugated polymers. 
First of all we note that our results for PAi and PPV 
are in good agreement with previously published ^ 85 * 92 ! 
GoWo results. In particular, for PAi we obtain a GW gap 
E g = 2.05 eV, to be compared with 2.1 eV (Ref. [55155] ) 
and 2.13 eV (Ref. |92J). These results are obtained for 
isolated chains of PA. In the case of a crystal, the gap 
shrinks" 84 to 1.8 eV due to interchain intera ctions. For the 
isolated chain of PPV, Rohlfing et aZPlSSl found E g = 3.3 
eV (Go Wo), which compares reasonably well with our 
G W result of 3.09 eV [see Tab. [TV]. In general, the 
overall shape (including the band widths) of the GW- 
corrected band structures for PA and PPV computed in 
this work is in excellent agreement with that of Ref. [85] . 

From a methodological point of view, the accuracy of 
G W corrections (based on the LDA electronic struc- 
ture) fo r organ ic molecules has been recently widely ad- 
dressedpSMEMSI] com p ar i n g with different implementa- 
tions of self-consistent GW and experimental results. 
Go Wo (LDA) is found to underestimate ionization po- 
tentials more than in extended systems, suggesting that 
a certain degree of self-consistency tends to improve on 
the results. Moreover, self-consistency is found to further 
lower the HOMO level and increase the fundamental gap. 
This leads to larger estimates of f3 max . 

In terms of a direct quantitative comparison with ex- 



periments, some issues have to be taken into account. 
First, electronic structure (and optical) measurements for 
polymers typically distinguish between crystalline grains 
and amorphous regions. Isolated chains are considered to 
resemble more (and to be used as rough models for) the 
amorphous regions. Clearly, a direct theory-experiment 
comparison may suffer from systematic errors (interchain 
interactions, medium polarization, electrostatic effects). 
These features generally tend to reduce the fundamental 
gap wrt that of the ideally isolated chain. With this in 
mind, we can compare data calculated here with experi- 
mental data from photoemission (PES) or scanning tun- 
neling (STS) spectroscopies. Rinaldi et al., measure 

d 97 

the electronic gap of PPV films (on a GaAs substrate) 
by means of STS. They were able to estimate E g ~ 3 
eV. Kemerik et alP^ used also STS and found the fun- 
damental gap of PPV [film deposited on Au(lll)] to be 
around 2.8 eV. All of these results are to be reasonably 
considered as lower bounds of the theoretical gap for the 
isolated PPV chain. 



C. Complex band structure: trends 

As can be directly inferred from the model described in 
App. [A] as well as from Ref. [5J, the important parame- 
ters that determine the behavior of /3(E) (and /3 max ) are 
the band gap E g , and the effective band widths of the 
states around the gap, given eg in terms of the hopping 
t\. While our model (see App. [A| includes a second pa- 
rameter t2 to describe the difference in the band widths of 
the frontier bands, considering that the ratio tijt\ ranges 
from 0.1 to 0.01 or less, corrections to the model due to 
t 2 are not particula rly releva nt for the cases studied here. 
According to Eqs. ( A4|A18 ), /3 max is mostly determined 
by the E g /t\ ratio. Even though this is just a simplified 
NN model, our numerical investigations suggest that the 
model is widely applicable (using a folding technique as 
in the case of PPI when needed) . 

In general the band gaps are expected to increase with 
the fraction of non-local exchange included in the hy- 
brid functional. For covalently bonded systems, the same 
trend is expected for the band widths. Numerically, in 
all our examples, while increasing the energy gap, HF 
increases also the band widths. The same trend is found 
for all the hybrid functionals we have investigated. Since 
both E g and t\ increase, it is not trivial to understand a 



priori which mechanisms would dominate. Indeed, it is 
very clear that the band gap opens more than the band 
width, leading to a clear trend of /3 max increasing when 
a larger fraction of exchange is included in the calcula- 
tion. The same trend is also found for the GW results, 
even though we had to extrapolate the CBS from the real 
band structure by using the model fitting (as described 
in details in Sec. IV A). 



In Fig. [9] we report a synthetic view of all the computed 
values of /3 max (times a/2, a being the polymer lattice pa- 
rameter), including all electronic structure methods for 
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FIG. 9: (Color online). Computed /3 max o/2 versus E g /ti, 
Eg being the band gap, t\ the effective hopping, and a the 
lattice parameter. All the polymers and XC functionals or 
methods are plotted. Black circles: PAi, squares: PA2; blue 
diamonds: PPV; orange triangles: PPL Open symbols refer 
to HF results. Dashed red line: /3 max according to the tight- 
binding model. 



PA (PAi, PA2 are shown; PAhf is not shown because 
almost ovcrimposed to the other PA geometries), PPV 
and PPL Those data are plotted against the E g /t\ value. 
The ideal curve from the tight-binding NN model is re- 
ported (dashed red line). The agreement between the 
computed and the modelled data is remarkable for all the 
cases studied. HF data are reported with empty symbols, 
showing in general a slightly worse agreement (as already 
discussed in Sec. 



IV C). On the basis of the above rela- 
tions, and according to the results reported in Fig. [9j we 
suggest the use of the model fitting to extrapolate infor- 
mation about /3 max from experiments able to investigate 
the electronic structure. This would allow for an indi- 
rect measure of the CBS and the related parameters (as 

/^max ) • 



D. Comparison with transport data 

As discussed in the introduction, in order to compute 
the f3 decay of the current (or conductance) in a metal- 
insulator-metal (MIM) junction, we need two different in- 
gredients: (i) the knowledge of the CBS for the infinitely 
long insulating region, and (ii) the position of the Fermi 
level of the MIM junction wrt the band structure of the 
insulator. This is depicted in Fig. [T] Assuming that at 
low bias the current is carried by the states close to the 
Fermi level of the junction, we basically need to know 
how these states decay into the insulating region (the 
polymer, in the present case). Once the CBS is known, 
either the Fermi level alignment is computed explicitly 
or it is estimated on the basis of physical considerations. 
While the direct calculation is feasible (but demanding) 
and in some cases necessary, it is also possible to give a 
first estimate of the Fermi level according to the s o-called 
MIGS (metal-induced gap states) theoryP^Hl As pro- 
posed by TersoffpZl if the metal DOS around the Fermi 
level is sufficiently featureless and the MIGS penetrate 



deep enough in the insulating region, the Fermi level 
of the metal-insulator junction is approximately pinned 
at the charge-neutrality level (often close to the midgap 
point) in order to avoid charge imbalance at the interface. 
The charge neutrality level can be easily identified from 
the CBS, as the energy where (3(E) reaches its maximum 
inside the gap. From this perspective, /3 max is a first esti- 
mate of the experimental (3 decay. In general, the Fermi 
level will move from the charge neutrality level, and the 
j3 decay will change accordingly. The physical reasons 
leading the Fermi level to shift are mostly related to the 
charge transfer (and dipole formation) at the interface. 
This explains why different chemical linking groups on 
the same molecule may lead to different values of (3. In 
general, the /3 max value computed from the CBS may be 
regarded as an upper bound for (3 and thus as a lower 
bound for the ability of the insulating layer to allow the 
current to tunnel through the junction. 

The above discussion stresses the fact that the ex- 
perimentally measured (3 does not in general depencP 
only on the electronic structure of the insulator (espe- 
cially the CBS), but also on the details of the interface 
(which determines the position of the Fermi level). For 
instance, recent measurements on alkanes (oligomers of 
PE) determined a (3 decay length of 0.71-0.76 A -1 for 
for -NH2 terminations J122H122I while it has been found 
in the range 0. 8-0.9 A - 1 (and more disperse) for thio- 
lated molecules ! 5 * 1 02 l EE] Recent calculations confirm this 
picture also for conjugated polymers connected to gold 
leads through different chemical groupsP In that work, 
the Authors have studied a number of oligomers includ- 
ing oligo-phenylenes (whose infinite polymer is poly-para- 
phenylene, PPP). In the case of the molecule connected 
to gold leads via thiol groups, after investigating the 
interface-DOS projected on the molecule, they find that 
the Fermi level aligns close to the midgap point (reason- 
ably the charge neutrality level considered here), off by 
few tenths of eV (^0.2-0.3 eV, the molecular gap being 
about 2 eV) at the LDA level. Since the basic unit of 
PPP (a phenyl ring) is the same as part of the monomers 
forming PPV and PPI, we assume the charge neutrality 
condition to be almost fulfilled if we were considering Au- 
PPV-Au and Au-PPI-Au junctions with thiol anchoring. 
We can thus estimate /3 max to be a good estimate of /?, 
keeping in mind that the a small deviation from midgap 
would slightly decrease f3. 

Recent experiments^!! reported j3 = 0.3 A -1 for 
oligomers of PPI connected to gold through thiols. 
This number is in very good agreement with the value 
/3 max = 0.29 A" 1 we found for PPI using PBE0 (see 
Tab. [y]). Acc ording to our findings for PA and PPV (see 
Tabs. III|IV ), GW should give results for (3 comparable 
to PBE0, but slightly larger. In order to compare with 
other existing (experimental and theoretica l) results for 
oligo-phenylenes (PPP in the infinite limit p | i03 | i 05ii06| 
we would need to address separately the issue of the 
Fermi level alignment (tending to decrease j3 wrt the ideal 
/3m ax ) and that of the phenyl twist-angle (which goes in 
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the direction of increasing 0), This will be the subject 
of future work. Coming t o the ca se of PA, we note that 
recent experimental results^ 7 * 108 * report /?- values of 0.22 
A -1 for molecules similar to oligo- acetylenes. While this 
number is in very good agreement with our GW (and 
PBEO) results for PAi (and consistent with PAhf, see 
Tab. |IIip , the large variability of (3 with the structural 
parameters of PA does not allow us to be conclusive on 
the assessment of the theory vs the experiment for this 
specific case. Nevertheless, our findings suggest that the 
use of PBEO or GW results, together with a proper de- 
termination of the Fermi level alignment, will provide a 
reasonable approximation. We also note that gap under- 
estimation and a priori assumption of the validity of the 
MIGS theory tend to partly cancel each other, and fortu- 
itous agreement of experimetally measured j3 values with 
LDA calculations may occur. 



VI. CONCLUSIONS 



while further developments along the lines of Ref. [70] 
are probably needed. 
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In this work we have computed from first principles the 
real and complex band-structure of prototype alkylic and 
conjugated polymer-chains using a number of theoreti- 
cal schemes, ranging from local and semilocal to hybrid- 
DFT, and GW corrections. The accuracy of these differ- 
ent methods has been evaluated and compared with ex- 
isting theoretical and experimental data, both in terms of 
the electronic structure and transport properties. From 
the CBS the (3 decay parameter, which governs non- 
resonant tunneling experiments through metal-insulator- 
metal junctions, can be computed. 

In doing so we have stressed the formal analogy of 
hybrid-DFT and GW (especially in the COHSEX formu- 
lation), and the interpretation of the hybrid-DFT elec- 
tronic structure as an approximation to the proper quasi- 
particle spectrum. We have also described in detail how 
to interpolate GW results by using a Wannier function 
scheme. In this case we have found that while the real 
band structure is always well interpolated, the CBS needs 
the self-energy real-space decay to be properly treated 
(off-diagonal representation and self-consistency of the 
wavefunctions) . 

We have numericaly investigated four polymers, 
namely poly-ethylene (PE), poly-acetylene (PA), poly- 
para-phenylcnc-vinylene (PPV), and poly-phenylcnc- 
imine (PPI). Our results compare well with the existing 
theoretical and experimental literature. Among the hy- 
brid functionals studied, PBEO results compare best with 
the Go Wo electronic structure. While the band gaps may 
still have a non-negligible deviation from GW, the agree- 
ment is remarkable on the CBS and (3 coefficient. The 
comparison with transport data (when available) is also 
very promising. This suggests PBEO as an efficient and 
reliable alternative to GW for these class systems, at least 
for transport properties. More generally, a systematic ap- 
plication of hybrid functionals to improve the accuracy 
of DFT-based electronic structure results is appealing, 



Appendix A: One-dimensional model 

In the case of one dimensional systems like conjugated 
polymers, numerical results for j3 can be rationalized in 
terms of a simple tight binding model as presented in 
Ref. [5]. In their work, Tomfohr and Sankey presented a 
two-band model which provides an analytical expression 
for the complex band structure (CBS) within the funda- 
mental energy gap. This model can also be used to fit 
the real band structure of realistic systems in order to 
evaluate the CBS analytically. The model is described 
in terms of two inequivalent sites (e a ,b) with nearest- 
neighbors (NN) hopping t\. These three parameters can 
be recast into E g (the fundamental gap), t\ and a further 
shift of the energy levels, which has no physical meaning. 
Moreover, it is shown that in this case /3 max (the maxi- 
mum value of /3(E) within the fundamental gap) depends 
only on E g /ti, according to: 



(3(E) a/2 



hi 



j(E) + v/7(£) 2 - 1 



(E-E V )(E C -E) 



= £ = 



2t\ 

E v + E c 



= 1 



1 (Ei 



(Al) 
(A2) 
(A3) 
(A4) 



All the details are given in Ref. [8] . We note however that 
this model is unable to reproduce any difference in the 
widths of the two bands (thus resulting in the CBS max- 
imum being located at mid-gap), while in our realistic 
simulations we typically find the LUMO bandwidth to be 
somewhat greater than that of the HOMO. We have then 
generalized the above to the following three-parameters 
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FIG. 10: Cartoon of the Hamiltonian adopted to fit the 



computed data, according to Eq. ( A5 \ 



model, which is more suitable to determine the physi- 
cal relevant quantities from our simulations. Extending 
the previous model, we add second NN interactions (with 
strength £2) between equivalent sites, t\ being the hop- 
ping between inequivalent sites as in the previous model. 
The model Hamiltonian is the following: 



H = 



R 



R 

[^a,R<R + ^1,R.+ A,R " 

R 

2 h K.R+A.R + ^6,R+lA,R 



■ cc + 



cc, (A5) 



where a, b indicate inequivalent sites and R is a cell index. 



The model is pictorially described in Fig. 10 The param- 
eters *i, * 2 are in general complex numbers. Taking t\, 
* 2 to be real for semplicity, the analytical expressions of 
the energy bands read: 



£1,2 (fc) = £ + 2i 2 x(k) ± [A 2 + 2t\ x(k) 
where, assuming e a > £(,, we have set: 



E v — ef, — 2* 2 ; E c — e a 
x(k) = 1 + cos(fca); 



2U 



E = \{E C 



E v ) ■ 
E„) 



E n 



(A6) 

(A7) 
(A8) 

(A9) 
(A10) 



In this picture E c _ v are the onsets of the valence and con- 
duction bands, a is the lattice parameter, k runs from — ~ 
to ~ . Under the condition |t 2 | *C |*i|, the band gap of 
the model is still at the Brillouin zone edge. A differ- 
enc choice of the relative phases of the model parameters 
would be needed to have the band gap at T (which is the 
case for PPI). 

It is also possible to express the results of the model 
in terms of more physical parameters such as the energy 
gap and the HOMO and LUMO band widths; 



W c = [A 2 +4< 2 ] 1/2 -A 



W v 



A 



M\ 1/2 



A 



4i 2 
4i 2 



(All) 
(A12) 



Because in the realistic calculations the bands of interest 
may cross other bands far from the k corresponding to 



the gap, we consider partial (W c , W v ) and not full band 
widths. These quantities are defined by the amplitude 
of the bands in a limited range of the BZ around the 
fundamental gap. This yields a much better agreement 
of the model bands with the calculated bands close to 
Eg, which is the energy range of interest. In order to 
extract the parameters E g ,t\, t 2 from our calculations we 
used the following relations (under the restriction that 
the band gap is direct and located at the Brillouin zone 
edge k = it /a): 



x = 
ti = 



t? = 



x(k ), 

4^- 

2^ 



- 2t 2 xq 



A) 2 - A ; 



(A13) 
(A14) 

(A15) 



Following Ref. [5], once we have parametrized the 
model, we can give an analytical expression for the CBS, 
which in our case reads: 



l(E) 



(E~E V ){E C -E) 
2(tj + 2Et 2 - 2i 2 E) 



+ 1, 



(A16) 



where Eq. ( Al ) connecting f) to 7 holds unchanged. The 



maximum of j(E) can be found analytically, and the re- 
sulting expression can be further simplified under the as- 
sumption |*2 1 <C |*i |; 



E r , 



7(£/w) ~ 1 + 



A 

*T 

Eg 
*1 



(A17) 



(A18) 



Appendix B: Computational details 

The CRYSTAL09 software package^ perform calcu- 
lations based on the expansion of the crystalline orbitals 
as a linear combination of a local basis set consisting of 
atom centered Gaussian orbitals. A 6-31G* contraction 
double valence (one s, two sp and one d shells) qual- 
ity basis sets have been selected to describe carbon and 
nitrogen atoms; the most diffuse sp (d) exponents are 
a c = 0.1687 (0.8) and a N = 0.2120(0.8) Bohr.- 2 The 
hydrogen atom basis set consists of a 31G* contraction 
(two s, one p shells): the most diffuse s and p exponents 
are 0.1613 and 1.1 Bohr. -2 The self consistent field pro- 
cedure was converged to a tole ranc e in the total energy of 
AE = 2 ■ 10~ 7 Ry per unit cellP^I Reciprocal space sam- 
pling was performed on Monkhorst-Pack grid with 12 k- 
points. The thresholds for the maximum and RMS forces 
(the maximum and the RMS atomic displacements) have 
been set to 0.00090 and 0.00060 Ry/Bohr (0.00180 and 
0.00120 Bohr). 

The calculations performed with Quantum- 
ESPRESSO (QE) adopt a plane waves basis set 
and norm-conserving pseudopotentials to describe the 
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ion-electron interaction. The kinetic energy cutoff 
has been set to 45 Ry for wavef unctions. For ionic 
relaxation, total energy LDA and GGA calculations use 
a Monkhorst-Pack grid of 8, 6, and 6 k-points for PE, 
PPV and PPI respectively (PA geometries are taken 
from the literature and not relaxed with QE). The 
convergence threshold on the atomic forces has been 
set to 10~ 3 Ry/Bohr. A minimum distance of 20 Bohr 
between chain replica is used. 

When performing GW calculations using SaX, the k- 
points grids have been made finer by using 50, 50 and 20 
k-points for PA, PE and PPV, respectively. The long- 
range divergence of exchange-like Cou lomb integrals is 
treated using a generaliz ed v ersiorP"^ of the approach 
given by Massidda et alP^ The same approach has 
also been used when performing hybrid-DFT calculations 
with QE. Note that other schemes to treat exchang e in 
one-dimensional systems have been proposedP^KIS] -phe 



Godby-Needs plasmon-pole modef^ has been used, set- 
ting the fitting energies at 0.0 and 2.0 Ry along the im- 
maginary axis, for all the cases. A kinetic energy cutoff 
of 6 Ry has been used to represent the polarizability P 
and the dynamic part of the screened Coulomb interac- 
tion W on a plane- wave basis, while a cutoff of 45 Ry 
has been used for the exchange operator. In order to 
converge the sums over empty states for the polarizabil- 
ity (self-energy), a total number of 288, 288, 288 (288, 
288, 608) states has been used for PA, PE, and PPV re- 
spectively. This corresponds to an equivalent transition- 
energy cutoffs of 53, 53, 35 eV (53, 53, 44 eV). Interchain 
distance has been increased to ~30 Bohr to control spu- 
rious interactions of periodic replica. The QP corrections 
are computed by evaluating the diagonal matrix elements 
of the self-energy operator (nk|E xc |nk), unless explicitly 
stated. 
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